/* Copyright (C) 2005, 2006, 2007, 2008, 2009, 2010 Mario Orsi
   This file is part of Brahms.
   Brahms is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by
   the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
   Brahms is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
   You should have received a copy of the GNU General Public License along with Brahms.  If not, see <http://www.gnu.org/licenses/>. */

// forceNbdStky.c -- functions for evaluation of the SSD "sticky" potential 

#include "definitions.h"

/* local (to this module) variable -- values from Fennel and Gezelter, J Chem Phys, 120, 9175 (2004) */
static const real v0 = 3.90 * cal_IN_J, // [kJ/mol]
  w0 = 0.07715;  
static const real halfv0 = 0.5 * 3.90 * cal_IN_J; // v0 / 2 [kJ/mol] 
static const real rL = 0.24, // [nm]
  rL1 = 0.275; // [nm]
const real rU  = 0.38, // [nm] 
  rU1 = 0.335; // [nm]

/*************************************************************************************************************************************
 * CalcStickyInteraction -- calculates SSD-SSD "sticky" pair interaction. Reference: Chandra & Ichiye, J Chem Phys, 11, 2701 (1999). *
 *************************************************************************************************************************************/
void CalcStickyInteraction( real *uij,           // out - potential energy
			    VecR *stickyForce,    // out - force
			    VecR *stickyTorqueij, // out - torque
			    VecR *stickyTorqueji, // out - torque
			    const VecR dr,        // in - distance
			    const RMat rMati, 
			    const RMat rMatj,
			    const real rr )       // in - r*r 
{
  real s, s1; /* switching factors */
  real r, sFactor, rUmr, CrUmrL;
  real rrr; /* cube( pair distance ) */
  real wij, wji, w;  /* tetrahedral sticky potential energy */
  real w1ij, w1ji, w1 = 0.; /* tetrahedral sticky potential correction energy */
  real fwijFactor, fwjiFactor, rbijzSq, rbjizSq, xxijMinusyyij, /* xij^2 - yij^2 */
    xxijMinusyyijDivrrr, /* (xij^2 - yij^2)/rrr */
    zijDivrr, /* z/rr */
    zijDivrrr, /* z/rrr  */
    zr6ij, /* zij/rbij - 0.6 */
    zr8ij, /* zij/rbij + 0.8 */
    zr62ij, /* (zij/rbij - 0.6)^2 */
    zr82ij; /* (zij/rbij + 0.8)^2 */
  real xxjiMinusyyji, /* xji^2 - yji^2 */
    xxjiMinusyyjiDivrrr, /* (xji^2 - yji^2)/rrr */
    zjiDivrr, /* z/rr */
    zjiDivrrr, /* z/rrr  */
    zr6ji, /* zji/rbji - 0.6 */
    zr8ji, /* zji/rbji + 0.8 */
    zr62ji, /* (zji/rbji - 0.6)^2 */
    zr82ji; /* (zji/rbji + 0.8)^2 */
  real ri, stickyScalingFactorij, stickyScalingFactorji;
  VecR drji, rbij, rbji, fibwij, fibw1ij, fibwji, fibw1ji, Fib_ij, Fib_ji, 
    FiSij, FiSji, Tibwij, Tibwji, Tibw1ij, Tibw1ji, Ti_ij, Tj_ji, switchGradient;
  
  VZero( fibw1ij );
  VZero( Tibw1ij );
  VZero( fibw1ji );
  VZero( Tibw1ji );
  // rbij is defined as the position of molecule j in the frame fixed on molecule i, so rbij = rMati * (-drij)
  VSCopy( drji, -1., dr );
  MVMul( rbij, rMati.u, drji ); /* rij = Ri * drji */
  //	    printf("*** rbij = % f % f % f\n", rbij.x, rbij.y, rbij.z);
  xxijMinusyyij = Sqr( rbij.x ) - Sqr( rbij.y );
  r = sqrt ( rr );
  rrr = r * rr;
  xxijMinusyyijDivrrr = xxijMinusyyij / rrr;
  zr6ij = rbij.z / r - 0.6;
  zr8ij = rbij.z / r + 0.8;
  zr62ij = Sqr( zr6ij );
  zr82ij = Sqr( zr8ij );
  zijDivrr = rbij.z / rr;
  zijDivrrr = zijDivrr / r;
  rbijzSq = Sqr( rbij.z );
  fwijFactor = - 6. * zijDivrr * xxijMinusyyijDivrrr;
  wij = 2. * xxijMinusyyijDivrrr * rbij.z; /* eq. (8) from Chandra & Ichiye, J Chem Phys, 11, 2701 (1999) */
  /* evaluate wij-part of sticky force on i */
  fibwij.x =   4. * rbij.x * zijDivrrr  + fwijFactor * rbij.x; /* (12) from Chandra & Ichiye, J Chem Phys, 11, 2701 (1999) */
  fibwij.y = - 4. * rbij.y * zijDivrrr  + fwijFactor * rbij.y; /* (13)[ditto for this and following numbers] */
  fibwij.z =   2. * xxijMinusyyijDivrrr + fwijFactor * rbij.z; /* (14) */
  /* evaluate wij-part of sticky torque on i */
  Tibwij.x =   4. * ( rbij.y * rbijzSq + rbij.y * xxijMinusyyij / 2. ); /* (21) */
  Tibwij.y =   4. * ( rbij.x * rbijzSq - rbij.x * xxijMinusyyij / 2. ); /* (22) */
  Tibwij.z = - 8. * VProd( rbij );                                      /* (23) */
  VScale( Tibwij, 1. / rrr );
  MVMul( rbji, rMatj.u, dr ); /* rbji contains dr_ji's coordinates in molecule-j's frame */
  xxjiMinusyyji = Sqr( rbji.x ) - Sqr( rbji.y );
  xxjiMinusyyjiDivrrr = xxjiMinusyyji / rrr;
  zr6ji = rbji.z / r - 0.6;
  zr8ji = rbji.z / r + 0.8;
  zr62ji = Sqr( zr6ji );
  zr82ji = Sqr( zr8ji );
  zjiDivrr = rbji.z / rr;
  zjiDivrrr =  zjiDivrr / r;
  rbjizSq = Sqr( rbji.z );
  fwjiFactor = - 6. * zjiDivrr * xxjiMinusyyjiDivrrr;
  wji = 2. * xxjiMinusyyjiDivrrr * rbji.z;
  w =  wij + wji ;
  /* evaluate wji-part of sticky force on i */
  fibwji.x =   4. * rbji.x * zjiDivrrr  + fwjiFactor * rbji.x; /* (12) */
  fibwji.y = - 4. * rbji.y * zjiDivrrr  + fwjiFactor * rbji.y; /* (13) */
  fibwji.z =   2. * xxjiMinusyyjiDivrrr + fwjiFactor * rbji.z; /* (14) */
  /* evaluate wji-part of sitcky torque on j */
  Tibwji.x =   4. * ( rbji.y * rbjizSq + rbji.y * xxjiMinusyyji / 2. );             /* (21) */
  Tibwji.y =   4. * ( rbji.x * rbjizSq - rbji.x * xxjiMinusyyji / 2. );             /* (22) */
  Tibwji.z = - 8. * VProd( rbji );                                                  /* (23) */
  VScale( Tibwji, 1. / rrr );
  if ( r > rL ) { /* switching */
    rUmr = rU - r;
    CrUmrL = Cube( rU - rL );
    s = Sqr( rUmr ) * ( rU + 2. * r - 3. * rL ) / CrUmrL;
    w *= s;
    VScale( fibwij, s );
    VScale( fibwji, s );
    VScale( Tibwij, s );
    VScale( Tibwji, s );
    sFactor = 6. * rUmr * ( rL - r ) / ( CrUmrL * r );
    VSCopy( switchGradient, sFactor, rbij );
    VVSAdd( fibwij, wij, switchGradient );
    VSCopy( switchGradient, sFactor, rbji );
    VVSAdd( fibwji, wji, switchGradient );
  }
  if ( r < rU1 ) { /* compute "w1"-part */
    /* w1ij */
    ri = 1. / r;
    w1ij = zr62ij * zr82ij - w0;
    fibw1ij.x = -2. * rbij.x * zijDivrrr;          /* (15) */
    fibw1ij.y = -2. * rbij.y * zijDivrrr;          /* (16) */
    fibw1ij.z =  2. * ( ri - zijDivrrr * rbij.z ); /* (17) */
    stickyScalingFactorij = zr62ij * zr8ij + zr6ij * zr82ij;
    VScale( fibw1ij, stickyScalingFactorij );
    Tibw1ij.x =   2. * rbij.y * ri;                /* (24) */
    Tibw1ij.y = - 2. * rbij.x * ri;                /* (25) */
    Tibw1ij.z = 0.;                                /* (26) */
    VScale( Tibw1ij, stickyScalingFactorij );
    /* w1ji */
    w1ji = zr62ji * zr82ji - w0;
    w1 = w1ij + w1ji;
    fibw1ji.x = -2. * rbji.x * zjiDivrrr;          /* (15) */
    fibw1ji.y = -2. * rbji.y * zjiDivrrr;          /* (16) */
    fibw1ji.z =  2. * ( ri - zjiDivrrr * rbji.z ); /* (17) */
    stickyScalingFactorji = zr62ji * zr8ji + zr6ji * zr82ji;
    VScale( fibw1ji, stickyScalingFactorji );
    Tibw1ji.x =   2. * rbji.y * ri;                /* (24) */
    Tibw1ji.y = - 2. * rbji.x * ri;                /* (25) */
    Tibw1ji.z = 0.;                                /* (26) */
    VScale( Tibw1ji, stickyScalingFactorji );
    if ( r > rL1 ) { /* switching, for rL1 < r < rU1 */
      rUmr = rU1 - r;
      CrUmrL = Cube( rU1 - rL1 );
      s1 = Sqr( rUmr ) * ( rU1 + 2. * r - 3. * rL1 ) / CrUmrL;
      w1 *= s1;
      VScale( fibw1ij, s1 );
      VScale( fibw1ji, s1 );
      VScale( Tibw1ij, s1 );
      VScale( Tibw1ji, s1 );
      sFactor = 6. * rUmr * ( rL1 - r ) / ( CrUmrL * r );
      VSCopy( switchGradient, sFactor, rbij );
      VVSAdd( fibw1ij, w1ij, switchGradient );
      VSCopy( switchGradient, sFactor, rbji );
      VVSAdd( fibw1ji, w1ji, switchGradient );
    }
  } /* end if ( r < ru1 ) */
  *uij = halfv0 * ( w + w1 );
  VAdd( Fib_ij, fibwij, fibw1ij ); /* i-frame */
  VAdd( Ti_ij,  Tibwij, Tibw1ij ); /* i-frame */
  VAdd( Fib_ji, fibwji, fibw1ji ); /* j-frame */
  VAdd( Tj_ji,  Tibwji, Tibw1ji ); /* j-frame */
  MVMulT( FiSij, rMati.u, Fib_ij );  /* FiSij = Ri^T * Fib_ij */
  MVMulT( FiSji, rMatj.u, Fib_ji );  /* FiSji = Rj^T * Fib_ji */
  VSub( *stickyForce, FiSij, FiSji );
  MVMulT( *stickyTorqueij, rMati.u, Ti_ij );  /* *stickyTorqueij = Ri^T * Ti_ij */
  MVMulT( *stickyTorqueji, rMatj.u, Tj_ji );  /* *stickyTorqueji = Rj^T * Tj_ji */
  VScale( *stickyForce,    halfv0 ); // 21/11/2006: fixed sign mistake
  VScale( *stickyTorqueij, halfv0 );
  VScale( *stickyTorqueji, halfv0 );
}
